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on-site approximation for the spin-orbit matrix elements. We have implemented the 
technique in the SIESTA[34] code, and show that it provides accurate results for the 
overall band structure and splittings of group IV and III-IV semiconductors as well as 
for 5d metals. 
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1. Introduction 

Computational methods in condensed matter physics are a powerful tool for predicting 
and explaining the most diverse properties of materials, nanostructures and small 
biological systems. Among an enormous plethora of methodologies, density functional 
theory (DFT) [TBI 125] has become a standard for simulations at the atomic and 
nanometric scale. Several numerical implementations of DFT are available and the 
details of the algorithm usually depends on the specific applications the method is 
designed for. These implementations can be categorized according to two different 
schemes. 

The first scheme divides DFT codes depending on the basis functions used, namely 
plane- waves or linear combination of atomic orbitals (LCAO). Plane waves are typically 
easier to implement and the convergence is determined by a single variational parameter, 
the energy cutoff. In contrast LCAO implementations are based on a tight binding 
description of the chemical bond. They are more cumbersome to implement and the 
variational principle is controlled by a set of parameters defining the basis functions. 
However these second methods are ideal for linear scaling since sparse Hamiltonian can 
be constructed by using strictly confined orbitals [T7] . 

The second classification takes into account whether the codes simulate both core 
and valence electrons or only the valence ones. In the first case the method is regarded 
as all electron since all the electronic degrees of freedom are treated on the same 
footing. This for instance allows one to perform fully relativistic calculations without 
any conceptual complication. All electron methods are remarkably accurate but have 
the drawback that the calculations are usually rather intensive and only relatively small 
systems can be tackled. In the second case the contribution of the core electrons is 
casted into pseudopotentials [27], which also can be constructed from DFT. This reduces 
drastically the number of electrons considered in the self-consistent simulation and much 
larger systems can be investigated. 

Spin-orbit interaction is a relativistic effect whose magnitude increases with the 
atomic number. Consequently it provides negligible contributions to the electronic 
structure of individual atoms and bulk materials made of light elements. However it has 
a significant impact over the physics of heavier elements, for instance 3d ferromagnetic 
materials. Spin-orbit can produce magnetic anisotropies of the order of 10 to 100 /xeV 
for bcc and fee Fe, Ni and Co [7J, therefore is the underlying mechanism for establishing 
the magnetic easy and hard axes and for controlling the shape of domain walls [6|. It 
is also the primary interaction responsible for most of the zero field splitting and other 
properties of magnetic molecules [T6j . 

In semiconductors spin-orbit interaction spin-splits the edges of the valence and 
the conduction band [10] and allows electrical manipulation of the spin-direction [10J. 
This last effect is of paramount importance for the growing field of spintronics [39J, 
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which has certainly added more impetus to the inclusion of the spin-orbit interaction 
in the description of the electronic structure. Spin-orbit coupling determines the 
spin-relaxation time of electrons in ordinary semiconductors and in semiconductor 
heterostructures, [21] and also plays an important role in the physics of diluted magnetic 
semiconductors. [35] Finally is worth mentioning that electron spin manipulation using 
spin-orbit interaction was recently demonstrated in the so-called spin-Hall effect [20] , a 
solid state version of the Stern-Gerlach measurement. 

It is therefore clear that spin-orbit interaction is becoming increasingly important 
for a number of applications, which also require the description of rather large systems. 
This calls for an efficient implementation of spin-orbit within pseudopotential LCAO 
based algorithms. Interestingly most of the mainstream LCAO codes such as SIESTA 
[34] . Onetep [33], Fireball [32] and Conquest [9] do not contain implementations of spin- 
orbit interaction in their present form. In contrast quantum chemistry packages such as 
Gaussian [12] or Turbomole [28] are not equipped for solid-state simulations. 

Therefore we have developed a general method for including spin-orbit interaction in 
conventional pseudopotential based LCAO DFT methods. This is not computationally 
demanding and hence it is fully adequate for large scale simulations. The method, 
although not suitable for highly accurate total energy calculations for which all electron 
plane-wave schemes cannot be matched, provides an excellent description of the effects 
of spin-orbit coupling on the electronic structure. Here we describe our implementation 
within the SIESTA program, [33] although the scheme is very general and it could be 
readily implemented in any LCAO pseudopotential codes with non-collinear electron 
spin functionalities. [291 CS] 

The paper is organized as follows: we first present the details of the method, our 
numerical implementation and numerical tests of one- and two-center integrals (section 
III). We then demonstrate the capability of the proposed scheme with predictions for 
group IV and III-V semiconductors and for 5d metals (sections IV and V, respectively). 

2. The On-site approximation 

2.1. Relativistic effects in pseudopotential methods 

Kleinman and Bylander have shown that the procedures for generation of non-relativistic 
pseudopotentials can be easily extended to account for relativistic effects. [231 M This 
relies on solving self-consistently the all-electrons Dirac equation for a single atom and 
in the extraction of a pseudopotentials Vj, where now j is the total angular momentum 
j — I ± |. The pseudopotential Hamiltonian can therefore be written as 



and includes all relativistic corrections to order a 2 , where a is the fine structure 
constant and \ j, mf) are total angular momentum states. This expression can be recast in 
a form suitable for non-relativistic pseudopotential theory by expressing \ j, mf) in terms 
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of a tensor product of the regular angular momentum states \l,m) and the eigenstates 
of the z component of the Pauli spin matrices [3 lj 



\j = l + -,mj) 



l,rrij + \) 



21+1 |">""j i 2 



1 v / V^Sf 1 K m i-|> 



-"" 7/ 1 y5pi^+i> 

Equation ([T]) can then be rewritten [I] as 

V = V sc + V 50 = 

= ^[^1 CT + ^ 5 °L-S] \l,m)(lM, (3) 

l,m 

where we use bold letters to indicate 2x2 matrices, with \ a representing the unit 
operator in spin space, 

L-S= 1 J l T z h h (4) 



2 \ L-L. — L? 



and, 



v, = 5 ^ T [(f+i)v; + . + ^..], 

It should be stressed that the scalar part V sc of the pseudopotential contains now not 
only the conventional non-relativistic pseudopotential, but also the scalar relativistic 
corrections. 

The vectors |Z, m), representing complex spherical harmonics, form a complete basis 
for the Hilbert space of the angular momentum operator L. It is a useful practice in 
solid state physics and quantum chemistry to replace them with real spherical harmonics 
\l,M), since the corresponding Hamiltonian is a real matrix. The change of basis is 
achieved by the following unitary transformation 

\l,M) =^= (\l,m) + (-ir\l,m)), 

\l,M) = -L (\l,m)-(-ir\l,m)), (6) 

which is valid for M > ( M = —M, m = —m). The case M = is simply given by 
\l,M = 0) = \l,m = 0). 

The pseudopotential operator V of equation ([31), is now written as 

V = V sc + V 5 ° = 

= ^[^U CT + ^ o r-S]|/,M)(/,M| . (7) 

LM 
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Finally the Kohn-Sham Hamiltonian [25] is expressed as a sum of kinetic energy T, 
scalar relativistic V sc , spin-orbit V 5 °, Hartree and exchange and correlation Y xc 
potentials 

H = T + V sc + V so + + V xc . (8) 
This Hamiltonian is therefore a 2 x 2 matrix in spin space 

fin fin 



H 



fin fin 



(9) 



whose non-diagonal blocks arise from the exchange and correlation potential whenever 
the system under consideration displays a non-collinear spin, and also from the spin- 
orbit potential. [29], [15] 

2.2. Spin-orbit in LCAO schemes: the on-site approximation 

LCAO methods expand the eigenstates \ip n ) of the non-collinear Kohn-Sham 
Hamiltonian over a set of localised orbitals \<f>i), 



n.i 

n 



i) 



(10) 



where % is a collective index for all the symbols required to describe uniquely a given 
orbital 

\4>i) = = \ R ^m) ® \ l i,Mi) . (11) 

Here both the radial and angular part of the wave function <f>i(r — dj) = (f \<pi) is 
centered at the position dj. Note that rij does not necessarily describe the principle 
quantum number only, but generally labels a set of radial functions associated to the 
angular momentum U according to the multiple zetas scheme. [30] 

The Kohn-Sham equation, H \ip n ) = E n \tp n ), is then projected onto such orbital 
basis set as 



1 1 



H ■ 



H 



n 



1 1 



till rp c 



71,7 

s 



0. 



(12) 



where Hff 



(falH^lh) and^ = (0 
and overlap matrix respectively. 

The spin-orbit term can then be calculated as 



) are the matrix elements of the Hamiltonian 



V 



= E<- 

k,l k ,M k 



(13) 



where index k indicates the atom on which the potential is centered, = V t (r — dk) 
and \lk, Mfc) is the spherical harmonic centered at the same atomic position dk- Equation 
([13]) involves a considerable number of three-center integrals. Inclusion of the spin- 
orbit interaction is therfore, in the LCAO approach, if straightforward, computationally 
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intensive. One possibility of reducing the computational effort is to construct fully 
non-local pseudopotentials. [21] We note however that the radial part of the spin- 
orbit pseudopotentials, is very short-ranged, resulting in matrix elements that 
decay quickly with the distance among the three centers. Thus we consider only 
matrix elements where the three components, both orbitals and the pseudopotential, 
reside on the same atom, and discard all the rest. This approximation simplifies the 
matrix elements of equation (1131) to one center integrals and drastically reduces the 
computational effort needed to account for spin-orbit effects. Then our approximated 
matrix elements are 



V 



so 



k,l h >0,M k 



b L 7 

L-t- — L? 



\l k ,M k ) (l^Mkll^Mj) 



(14) 



since the L a operators leave each Zj subspace invariant. The angular part of these on-site 
matrix elements can be calculated analyticallj||| 



(l,2\Lzp\l,Mj) 
{l,2\Lf\l,Mj) 
(l,3\L T \l,Mj) 
(1,311^11, Mj) 



± 



(I, Mi\L z \l, Mj) = - iMi5(Mi + Mj = 0) 



± 



± 
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£ This formula was not correctly written in the previous version of the manuscript. We wish to thank 
Hyungjun Lee from Yonsci University for kindly drawing our attention to this point. The correct 
formulae were in any case used in the code from the very beginning and, hence all the results of the 
simulations presented in the article are correct and remain the same. 
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Some useful symmetries of the matrix elements of the Hamiltonian and of its spin- 
orbit part should also be highlighted. Since the Hamiltonian is hermitian the matrix 
elements satisfy the relation 

Elf = (Ejf)*. (16) 

Moreover it is also easy to show that all the terms in the Hamiltonian satisfy a spin box 
hermiticity, i.e. 

Elf = (Elf)*, (17) 
except for the spin-orbit contribution which is spin box anti-hermitian 

Elf = -(Elf)*. (18) 
This property has important consequences for the calculation of the total energy. 

2.3. Density matrix and total energy 

The charge density can also be written in terms of the LCAO basis as 

n 

= J2 Ur~ cj>*(f- d s ) Pij (19) 
hi 

where /„ represents the occupation of the Kohn-Sham eigenstate ip n (f) and is a 2 x 2 
matrix containing the products of wave-function coefficients, whose components are 

Pif = J2f-<i c i*r ( 2 °) 

n 

The electronic contribution to the total energy may be expressed as a sum of a band- 
structure (BS) contribution plus double- counting corrections. The BS contribution can 
be written in the LCAO basis as 

E e S = ^fn Mm n ) = Jf, Hifp^i ( 21 ) 
n i,j,a,a' 

which may also be expressed as 



+2 Re [(V^'-Vf 3 ") (^r]} , (22) 

where we have isolated the non-diagonal contributions in spin space. These arise only 
from the spin-orbit interaction and from the exchange and correlation potential whenever 
spin non-collinearity is present. The spin-orbit contribution to the total energy therefore 
is 

,.;<o Tr^V-V„. (23) 
since there are no double- counting terms. 
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2.4- Forces and stresses 

One important consequence of the on-site approximation is that spin-orbit does not 
give rise to an explicit contribution to forces and stresses, even though an implicit 
contribution due to the modification of the self-consistent wave-function is always 
present. According to Hellmann-Feynman theorem, [13] the spin-orbit contribution 
to the forces excerpted upon an atom centered at position d\. is obtained by simply 
differentiating the energy with respect to the atomic coordinates of the atom, 

- F™=V k E so = 

hi 

where both the i and j orbitals are centered at the same atomic position dk and 
Vfe = V^ fe . However both contributions to the spin-orbit forces in equation (|24|) vanish. 
The first one is identically zero since the one-center integrals do not depend on the 
atomic position, 

V fe ( J R n ^|^f|it: n . > ,.)=0. (25) 
The second term may be rewritten as 

ij 

where £^°' CTfJ are the components of the spin-orbit contribution to the energy-density 
matrix 

e S° = \Y. O^ 1 V 2? Pmd + Pa V L° S-)). (27) 

However, since 8 is antisymmetric with respect to the orbital indices, in contrast to 
the overlap matrix that is symmetric, the second term vanishes as well. 

In a similar way one can demonstrate that the spin-orbit interaction in the on-site 
approximation does not introduce any additional contribution to the stress. 



(24) 



3. Numerical tests on one- and two-center integrals 

The validity of the on-site approximation relies on the fact that two and three center 
integrals are much smaller than the one-center ones, which are kept in the simulation. 
Among those two and three center matrix elements the two-center integrals 

V%(R) = (MR) \Vi-° (0) HO) ■ 5(0) |0,(O)), (28) 
are expected to have the largest absolute value. An excellent test consists of calculating 
these matrix elements for a given material along the direction that joins one atom with 
its nearest neighbours, as a function of the distance R between the two centers. Then, 
R = provides the on-site matrix elements, while if R equals the nearest neighbours 
distance, the calculation describes the desired two-center integrals. We have performed 
such test for a representative semiconductor, Si, and a representative 5c? metal, Pt. 
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Figure 1. Two-center integrals V^_ 2 .(i?) and V*_ x (i2) (solid and dashed lines, 
respectively) for silicon along the (111) direction, as a function of the distance R 
between the two centers. The arrow indicates the distance of the nearest-neighbour 
atoms. 



For silicon, the valence electrons include s- and p-orbitals only and therefore we 
consider the following matrix elements 

V^ X (R) = (cj> Py (R) \V p so (0)L z (0) 1^(0)) , 

V Z Z _ X (R) = (<f> Pz (R) \V p so (0) L,(0) |0 fc (O)) . (29) 

For R =0, the first of these matrix elements reduces to the on-site integral Vy_ x (0), 
while the second is zero by symmetry. Fig. [T] shows the matrix elements as a function of 
R along the direction (111). The matrix elements, evaluated at the nearest neighbour 
distance are considerably smaller than the on-site integral, with the V(i?)/V(0) ratio 
being ~0.03 and ~0.08 respectively for y — x and z — x. 

In the case of Pt we have considered not only the 5d but also the Qp orbitals. 
Platinum crystalizes with an fee structure and a nearest neighbours distance of 5.4 a.u. 
We compute the same p-matrix elements as for silicon, as well as the two following d 
integrals 

V'aW-vuW = (^-v< R ) \V d SO (0)L z (0) |^(0)> , 

V^_ M (i2) = (<h*-r*(R) \V d so (0) L z (0) |0^(O)) . (30) 

Fig. [2] shows these matrix elements as a function of R. For the rf-type matrix 
elements (figure EK) the decay with the distance between the centers is rather dramatic 
and the two-center matrix elements are about 10 -4 times smaller than the on-site 
integrals. The case of the p-integrals (figure[2b) is similar to that of Si with a V(i?) /V(0) 
ratio of ~ 0.08. For the specific case of Pt however these p-integrals are expected to 
contribute little to any physical quantities since they correspond to unoccupied states. 

Same tests for other materials give similar results. We therefore conclude that the 
on-site approximation provides accurate results for heavy transition metals and good 
quantitative estimates for semiconductors. For this latter we estimate an error due to 
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Figure 2. (a) Two-center integral Vj£a_ 3 _ V z z2 _ xz (R) (solid and dashed lines, 

respectively) and (b) V Z _ X (R) and VJ_ x (ii) (also, solid and dashed lines, respectively) 
for platinum along the (110) direction as a function of the distance R between the two 
centers. The arrow indicates the distance of the nearest-neighbour atoms. 



the neglecting of high-center integrals not larger than 10%, and on average of the order 
of 3-6%. 

4. Spin-orbit in group IV and III-V semiconductors 

We have thus demonstrated that in general our on-site approximation simplifies the 
inclusion of spin-orbit effects in LCAO DFT codes. We have thus implemented the 
method in SIESTA, a LCAO code able to simulate non-collinear arrangements of 
spins, [2H US] that in addition reads relativistically generated pseudopotentials in the 
form required by equation ([TJ). 

The numerical implementation is rather simple since the spin-orbit contribution 
to the Hamiltonian does not depend on the electron charge density and therefore does 
not need to be updated during the self-consistent procedure. This drastically reduces 
the computational overheads, which are almost identical to those of a standard non- 
collinear spin-polarised calculation. Here we present a series of test cases for the 
band-structures of group IV and III-V semiconductors, obtained with the local spin 
density approximation and norm- conserving relativistic pseudopotentials. Special care 
was taken in the generation of the pseudopotentials and of the basis sets and in the 
choice of the parameters that control the numerical accuracy of real and reciprocal 
space integrals. [34] In particular the basis set was highly optimised following the scheme 
proposed in references [ [T9l |3]], from which we have borrow our notation for multiple 
zeta basis sets SZ, SZP, DZ, DZP, TZ, TZP, TZDP, TZTP, TZTPF. 

The introduction of the spin-orbit interaction lifts specific degeneracies in the band- 
structure of a material. In particular, for diamond and zincblende semiconductors the 
six fold degenerate valence band r^ 5 , at T splits in two. The first is four-fold degenerate 
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Figure 3. Band structure of GaAs calculated with the on-site approximation to the 
spin-orbit interaction. 

Tg (the heavy and light hole bands), and the second is only two fold degenerate Tj 
(the spin-split-off band). This energy splitting Ao (T^ 5 — > T^,T^), is the hallmark of 
the effects of spin-orbit interaction in the band-structure of these materials. Other 
commonly measured energy splittings are called A' (Tl 5 — > r^Ty), Ai (Lg — > L v Ab) L§) 
and A; (L§->2£ b ,L§). 

We show in Fig. [3] the band-structure of the canonical III-V semiconductor, 
GaAs. The figure also defines graphically the energy splittings described in the previous 
paragraph. We note first that the band-structure closely resembles that obtained with 
other methodologies, and also agrees rather well with the experimental data. |38| except 
for the characteristic LDA underestimation of the semiconductor gap, which is further 
reduced because of the spin-orbit energy splitting. Therefore, in all the tests that follow, 
we have avoided narrow gap semiconductors, which are usually predicted to be metals 
by LDA. [36] 

4-1. Band structure of Si 

The spin-orbit energy splittings of silicon at the high symmetry points of the diamond 
lattice are presented in Table [1] and Fig. H] for increasingly complete basis sets. 

Our results are in extremely good agreement with the theoretical and experimental 
data available in the literature. [38] We note that a DZP basis set, which is usually 
assumed to be the minimal basis needed to obtain reasonably converged results, already 
provides accurate predictions for the energy splittings. The Aq split is somehow an 
anomaly, and in general we find that the splittings of the conduction bands are not as 
well described as those of the valence bands. 
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Basis 


E T feV) 


An (eV) 

1 — MJ V v y 




Ai feV) 

i — *i ^ v y 


Ai (eV) 


SZ 


-214.8 


0.051 


0.025 






DZ 


-215.0 


0.068 


0.157 


0.027 


0.105 


TZ 


-215.2 


0.068 


0.002 


0.023 


0.078 


SZP 


-215.8 


0.042 


0.787 


0.024 


0.032 


DZP 


-216.0 


0.044 


0.647 


0.024 


0.047 


TZP 


-216.0 


0.045 


0.696 


0.025 


0.050 


TZDP 


-216.0 


0.045 


0.604 


0.026 


0.043 


TZTP 


-216.0 


0.045 


0.593 


0.026 


0.044 


TZTPF 


-216.1 


0.044 


0.615 


0.025 


0.030 


REF 




0.044 


0.04 


0.02 


0.03 



Table 1. Spin-orbit energy splittings of bulk Si calculated for increasingly complete 
basis sets. Et is the total energy A , A' , Ai and are the splittings as denned in 
the text. REF corresponds to the reference values, experimental whenever possible, as 
described in the literature. [55] 




■215 -215.5 -216 

Total Energy (eV) 



Figure 4. Convergence of the spin-orbit energy splittings of bulk Si with the size of 
the basis set. From top to bottom panels, we present respectively A , Ag, Ai and A[ 
as functions of the total energy associated with each set. The points correspond to SZ, 
DZ, TZ, SZP, DZP, TZP, TZDP, TZTP and TZTPF, respectively. 

We note that Kohn-Sham eigenvalues should not be associated to single particle 
excitation energies, since the former are simply the Lagrangian multipliers leading to 
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Basis 


a (A) 


BM (GPa) 


E r feV) 


DZP 


5.392 


98.2 


5.480 


TZP 


5.389 


98.8 


5.485 


TZTP 


5.388 


98.2 


5.491 


PW 


5.384 


95.9 


5.369 


LAPW 


5.41 


96 


5.28 


EXP 


5.43 


98.8 


4.64 



Table 2. Structural parameters of bulk Si for several basis sets, a is the lattice 
parameter, BM the bulk modules and E c the cohesive energy. PW refers to a 50 
Ryd cutoff plane- wave calculation. [34j LAPW corresponds to an all-electrons linear- 
augmented-plane-wave calculation [14] and EXP to the experimental values. [22] 

the Kohn-Sham equations. This is valid for both the valence and the conduction 
bands. However it is commonly accepted that DFT band structures are a good first 
approximation to single particle energies of occupied states. The conduction bands 
are somehow different since such states do not contribute to either the total energy 
and the density matrix, and therefore are irrelevant in DFT. For this reason a stark 
disagreement in the conduction band splitting should not be surprising. Moreover, the 
systematic overestimation of the A' Q splitting is related to the underestimation of the 
bandgap. This produces an erroneous enhancement of the hybridisation between the 
orbitals forming the conduction bands with those forming the valence one, with a net 
overestimation of the spin-orbit splitting. It is therefore expected that corrective schemes 
to the bandgap may also correct the spin-orbit splitting of the conduction bands. 

It is also important to note that basis where the radial component varies sharply 
around the origin should be avoided. These in fact are difficult to integrate in the range 
where the spin-orbit pseudopotential is appreciable and brings considerable numerical 
instability to the evaluation of the matrix elements. 

Finally we have calculated the Si bulk parameters for different basis sets in order to 
check that the inclusion of spin-orbit interaction does not change significantly the LDA 
results. A summary of all the computed structural parameters is presented in Table |2] 

4 ■ 2. Ill- V semiconductors 

We have further tested our method by calculating the various energy splittings of several 
group IV and III-V semiconductors such as Ge, GaAs, AlAs and AlSb, i.e. of those with 
a reasonably large bandgap. For all of them we have found again that a DZP basis set 
provides essentially converged results. This is demonstrated in table [3] where we show 
that the splittings calculated with a DZP basis agree rather well with other theoretical 
estimates and with experimental values. Also in this case Aq is the exception for the 
same reasons explained before. 
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5. 5c? metals: Au and Pt 

Since the spin-orbit interaction is a relativistic effect, it is expected to increase with 
the atomic number. Therefore, metals from the fifth row of the periodic table are an 
ideal test ground for our method. Among those, gold and platinum are specially good 
candidates, since the first is a closed-d shell noble metal while the d-bands of the second 
have considerable weight at the Fermi energy. Moreover, a number of ab initio spin- 
orbit calculations are available [8] which demonstrate that the spin-orbit interaction is 
essential for the correct description of their band structures. To simulate these two 
materials, we have again used the LDA approximation for exchange and correlation 
potential, and constructed an optimised set consisting of two atomic wave functions in 
each of the s-, p- and d-channels. 

In Figs. [5] and [6] we present the band-structures of Pt and Au, calculated at 
the theoretical lattice constant, with and without taking into account the spin-orbit 
coupling, the figures show that s-bands are not modified by spin-orbit, while p- and 
d-bands suffer strong modifications, specially whenever two bands cross each other. 
Moreover, we find that the spin-orbit interaction lifts degeneracies at the band crossings 
as expected. 

We summarise in Table H] the bulk lattice constant and the energy of some selected 
bands at the T point, that we define graphically in Fig. [5j These energies are in good 
agreement with other much more expensive methods, like Plane-wave pseudopotential 
(PWSCF)[8] and relativistic full-potential Korringa-Kohn-Rostoker (KKR)[37] methods 
or augmented plane-wave (APW) approaches for solving Dirac equation. [TT] We note 
that the differences range in the order of only a few per cent. 



Material 


A (eV) 


A' (eV) 


Ai (eV) 


Ai (eV) 


Ge 


0.2959 


0.3783 


0.1545 


0.356 


REF01 


0.287 


0.200 


0.184 


0.266 


GaAs 


0.3573 


0.3006 


0.1857 


0.319 


REF[38] 


0.34 


0.26 


0.23 


0.11 


AlAs 


0.3073 


0.0762 


0.1636 


0.118 


REF[26] 


0.28 


0.04 


0.20 




AlSb 


0.6847 


0.1752 


0.3440 


0.307 


REF[38] 


0.75 


0.1 


0.4 


0.09 



Table 3. Spin splittings for several III-V semiconductors as calculated with a DZP 
basis set. REF correspond to reference values, experimental whenever possible, as 
described in the literature. 4. 38, 26 
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Figure 5. Band-structure of platinum obtained at the theoretical equilibrium lattice 
constant. The dashed line is for a spin-orbit calculation, while the solid line is obtained 
when the spin-orbit coupling is not included. The figure also provides a graphical 
definition of the energies at the F point, £i, of the bands that are closest to the Fermi 
energy. 



6. Conclusion 

We have presented a simple and effective method for including the spin-orbit interaction 
in standard LCAO DFT calculations, which is based on relativistic norm-conserving 
pseudopotentials and on an on-site approximation to the spin-orbit matrix elements. 
The method is computational undemanding and extremely simple to implement in 
standard LCAO DFT codes, such as SIESTA. Importantly the on-site approximation 
does not introduce additional contributions to both forces and stress. We have then 
presented a series of tests for group IV and III-V semiconductors and for 5d metals. 
The overall structural parameters do not change with respect to standard non-relativistic 
LDA calculations, and are in general good agreement with reference data. The spin- 
orbit splittings of the band structures are also in good agreement with those obtained 
with more computationally intensive DFT methods. 

The good results obtained for the electronic structures and structural parameters 
make our method very attractive for large scale simulations where spin-orbit coupling is 
relevant. This is for instance the case of semiconductors heterostructures and quantum- 
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Figure 6. Band structure of gold obtained at the equilibrium theoretical lattice 
constant. The dashed line is for a spin-orbit calculation, while the solid line is obtained 
when the spin-orbit coupling is not included. The figure also provides a graphical 
definition of the energies at the T point, of the bands that are closest to the Fermi 
energy. 



transport calculations [2] where spin- mixing effects are important. 
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